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Abstract 

Standard methods for describing the intensity distribution of mechanical and acoustic wave 
fields in the high frequency asymptotic limit are often based on flow transport equations. Common 
techniques are statistical energy analysis, employed mostly in the context of vibro-acoustics, and ray 
tracing, a popular tool in architectural acoustics. Dynamical energy analysis makes it possible to 
interpolate between standard statistical energy analysis and full ray tracing, containing both of these 
methods as limiting cases. In this work a version of dynamical energy analysis based on a Chebyshev 
basis expansion of the Perron-Frobenius operator governing the ray dynamics is introduced. It is 
shown that the technique can efficiently deal with multi-component systems overcoming typical 
geometrical limitations present in statistical energy analysis. Results are compared with state-of- 
the-art hp- adaptive discontinuous Galerkin finite element simulations. 

1 INTRODUCTION 



Predicting the wave energy distribution of the vibro-acoustic response of a complex mechanical system 
to periodic excitation is a challenging task, especially in the mid-to-high frequency regime. Standard 
numerical tools such as finite element methods become inefficient, and ray or thermodynamic approaches 
^ 'are often employed to model the wave energy flow through the structure. Popular methods are Statistical 
Energy Analysis (SEA) [TJ [2j [3], in which the mean energy flow between subsystems is assumed to be 
proportional to the energy gradient, and the ray tracing technique, in which the wave intensity distribution 
is determined by summing over contributions of a potentially large number of ray paths [HI SJ |5]. 

Ray tracing and SEA both predict mean values of the energy distribution and omit information about 
wave effects such as interference or diffraction. Both methods are therefore expected to hold in the high 
frequency (or small wavelength) limit. SEA is in fact a low resolution ray tracing method [SI [7] leading to 
small numerical models compared to ray tracing. This efficiency saving comes at a price, however: SEA 
has no spatial resolution of the energy distribution within subsystems and becomes unreliable whenever 
long range correlations in the ray dynamics are present. The recently developed Dynamical Energy 
Analysis (DEA) [7] provides a tool which interpolates between SEA and a full ray tracing analysis and can 
overcome some of the problems mentioned above at a relatively small computational overhead. DEA thus 
enhances the range of applicability of standard SEA and gives bounds on the range of applicability of SEA. 
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Related methods have been discussed previously in the context of wave chaos [9] and structural dynamics 
|10j . In particular Langley's Wave Intensity Analysis (WIA) [Tfl IT2] and Le Bot's thermodynamical 
high frequency boundary element method p3l E] include details of the underlying ray dynamics. The 
approach employed here differs from these methods by considering multiple reflections in terms of linear 
operators. Representing these operators in terms of basis function expansions then leads to SEA-type 
equations. 

In this work a reformulation of DEA is presented, which is based on a Chebyshev basis function 
representation; this leads to considerable improvements compared to previous attempts based on an ex- 
pansion in terms of a Fourier basis [7]. Both Chebyshev and Fourier basis expansions of smooth functions 
share similar exponential convergence properties [15J. The main advantages of using a representation 
in terms of Chebyshev polynomials include that the requirement for periodic boundary conditions can 
be dropped, allowing for much more freedom in the choice of approximation regions. In addition, a 
Chebyshev expansion gives rise to more efficient quadrature rules for numerically calculating the arising 
integrals when constructing the linear operators considered in DEA. Using a Chebyshev basis leads to a 
natural choice of quadrature, namely Gauss-Chebyshev, which is optimal for polynomial-type integrands 
and naturally incorporates the orthogonality weighting term in the Chebyshev basis function represen- 
tation. In order to take full advantage of this feature it is necessary to formulate the problem in terms of 
the final position and momentum of a given ray, and map back to its initial point. This is in contrast to 
previous work on DEA by Ref. [7J, where the rays were defined by their endpoints. The strengths of the 
newly reformulated DEA are evident in the applications considered, where due to improved efficiency it 
has been possible to model mult i- component systems with variable wave-numbers for the first time. The 
method is verified numerically by comparing DEA results with state-of-the-art finite element software 
for a range of parameter values. 

The remainder of the paper is structured as follows. In Section [21 the ray tracing approximation is 
discussed and related to the Green function using short wavelength asymptotics. In Section [3j the concept 
of phase-space operators is introduced and their representation in terms of boundary basis functions is 
discussed. In Section H] the implementation of DEA is detailed along with its links with SEA. The 
finite element formulation used for verification of the results is also briefly discussed. In Section [5j a 
variety of coupled two-cavity configurations are discussed and the results compared against finite element 
computations. Finally larger multi-cavity configurations are considered. 



2 WAVE ENERGY FLOW IN TERMS OF THE GREEN 
FUNCTION 



It is assumed that the system as a whole is characterized by a linear wave equation describing the overall 
wave dynamics including damping and radiation in a finite domain Q C M d , d = 2 or 3. In this work 
only stationary problems with continuous, monochromatic energy sources are considered. We split the 
system into Nq subsystems and consider the scalar wave equation for acoustic pressure waves in each 
homogeneous sub-domain with local wave velocity q, i — 1, ...,Nq and f2 = ILJi Extensions to 
more complicated systems with different wave operators in different parts of the system can be treated 
with the same techniques as long as the underlying wave equations are linear, see the discussion in Ref. 



The general problem of determining the response of a system to external forcing with angular fre- 
quency w at a source point r G Qq can then be reduced to solving 




(1) 
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with H = —A, G represents the Green function, r 6 f2j is the solution point and S is the Dirac delta 
distribution. Furthermore, ki = cu/ci + ifii/2 is a complex valued wavenumber, where the imaginary part 
represents a subsystem dependent damping coefficient /ij. Throughout this work we take i = a/— 1 unless 
used as a subscript, in which case it is an index over the number of subsystems. The wave energy density 
induced by the source is then given as 

\G(r,r ;u)\ 2 

e(r, r ; u) = - 2 , (2) 

QiCj 

for r G Qi where Qi is the density of the medium in The linear wave operator H can naturally be 
associated with the underlying ray dynamics via the Eikonal approximation; for more detailed derivation, 
see Ref. [TBI EH]- Using small wavelength asymptotics, the Green function in equation ([[]) may be 
written as a sum over all classical rays from r to r for fixed kinetic energy of the hypothetical ray 
particle. One obtains [HI [17] 

G(r, ro; W ) = l +1)/2 £ Aj^Wn (3) 

j:ro— >r 

where Rj is the length of the ray trajectory between r$ and r including possible reflections on boundaries. 
The amplitudes A, may be written as a product of three terms as in Ref. [7] due to damping, mode 
conversion and reflection/transmission coefficients, and geometrical factors. The phase index Vj contains 
contributions from the reflection/transmission coefficients at interfaces between subsystems and from 
caustics along the ray path. 

Analogous representations to ([3]) have been considered in detail in quantum mechanics [18] and are 
also valid for general wave equations in elasticity, see Ref. [5] for an overview. In the latter case G 
becomes matrix valued. Note that the summation in equation is typically over infinitely many 
terms, where the number of contributing rays increases (in general) exponentially with the length of the 
trajectories included. This gives rise to convergence issues, especially in the case of low or no damping 

M, 

The wave energy density ([2]) can now be expressed as a double sum over classical trajectories and 
hence 

e(r,r ;oo) =C £ ^^^[r^r^W.-u^^ 

j,j':ra->r (4) 

= C[p{r, r ; u) + off-diagonal terms], 

with C = 7r 2 / (gicj(27iY d+1 ^) . The dominant contributions to the double sum arise from terms in which 
the phases cancel exactly; one thus splits the calculation into a diagonal part 

p(r,r ;w)= £ \Aj\ 2 (5) 

j:ro—¥r 

where j = j', and an off-diagonal part. The diagonal contribution gives a smooth background signal 
and the off-diagonal terms give rise to fluctuations on the scale of the wavelength. The phases related 
to different trajectories are (largely) uncorrelated and the resulting net contributions to the off-diagonal 
part are in general small compared to the smooth part, especially when averaging over frequency intervals 
of a few wavenumbers. 

It has been shown in Ref. [7J that calculating the smooth diagonal part ([5]) is equivalent to a ray 
tracing treatment. That is, the smooth part of the energy density can be described in terms of the flow 
of fictitious non-interacting particles emerging from the source point r$ uniformly in all directions and 
propagating along ray trajectories. This makes it possible to relate wave energy transport with classical 
flow equations and thus thermodynamical concepts, which are at the heart of an SEA treatment. In DEA 
the classical flow is expressed in terms of linear phase space operators as detailed in the next section. 
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3 LINEAR PHASE SPACE OPERATORS AND DEA 

3.1 Phase space operators and boundary maps 

A brief outline of the derivation of the DEA flow equations is now given, for details see Ref . [7] . We adopt 
a purely kinetic viewpoint based on the interpretation that rays are trajectories of particles following 
Hamiltonian dynamics as detailed in Section 2 of Ref. [16]. Here the time dependence of a density of 
ray trajectories (or particles) p is known to satisfy the Liouville equation 

|W) + ^-Vx(p(X,r)) = 0, (6) 

where X = (r, p) denotes the phase space coordinate with position r and momentum p. The propagator 
for the Liouville equation is the linear phase space operator C T (X,Y) = 5{X — (p T (Y)), known as a 
Perron- Frobenius operator in dynamical systems theory [19], and hence we may write 

p(A,r)= / C T {X,Y)p Q {Y)dY. (7) 
Jp 

Here the phase space flow f T (Y) gives the position of the particle after time r starting at Y — (r',p') 
when r = 0. Furthermore, po denotes the initial ray density at time r = and the domain of integration 
is over the whole of phase space P = Q x M. d , where the integration over M. d takes care of the momentum 
coordinates p. 

Consider a source localized at a point r emitting waves continuously at a fixed angular frequency 
u>. Standard ray tracing techniques estimate the wave energy at a receiver point r by determining the 
density of rays starting at r and reaching r after some unspecified time. This may be written in the 
form 

poo r p 

p(r,r ,uj)= / w(Y,T)C T (X,Y)p (Y,uj)dYdpdr, (8) 



with initial density po{Y,u) = 5(r' — ro)S(kQ — H(Y)), where H = \p\ 2 is the Hamilton function for the 
wave operator H and fco is the wave number at the source point as defined in Eqn. ([1]). It can be shown 
that equation (jSJ) is equivalent to the diagonal approximation (jSJ) |7J. A weight function w is included 
to incorporate damping and reflection/transmission coefficients. It is assumed that w is multiplicative, 
(w(ip T1 (X), t 2 )w(X, fx) = w(X, Ti + r 2 )), which holds for standard absorbtion mechanism and reflection 
processes [T7] . 

In order to solve the stationary flow problem flH]) a boundary mapping technique is employed. For 
the time being let us consider a problem with a single (sub-)system Q = Q\ with boundary T. The 
boundary mapping procedure involves first mapping the ray density emanating continuously from the 
source onto the boundary T. The resulting boundary layer density p^ is equivalent to a source density 
on the boundary producing the same ray field in the interior as the original source field after one 
reflection. Secondly, densities on the boundary are mapped back onto the boundary by a boundary 
operator B(X s ,Y s ;u) = w(Y s )5(X s — where X s = (s,p s ) represents the coordinates on the 

boundary. That is, s parameterizes V and p s G B d ~^ 1 denotes the momentum component tangential to 

r at s for fixed H(X) = \p\ 2 , where -B^J 1 is an open ball in of radius \p\ and centre s. Likewise, 
Y s = (s',p' s ) and <fi^ is the invertible boundary map. Note that convexity is assumed to ensure W is well 
defined; non-convex regions can be handled by introducing a cut-off function in the shadow zone as in 
Ref. [13] or by subdividing the regions further. 
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The stationary density on the boundary induced by the initial boundary distribution p r j (X s , u) can 
then be obtained using 

oo 

PrM = ^£>)pp 0) (u;) = (I-BM)- 1 ^), (9) 

n=0 

where B n contains trajectories undergoing n reflections at the boundary. The resulting density distri- 
bution on the boundary p-p(X s ,u) can then be mapped back into the interior region. One obtains the 
density (jSJ) after projecting down onto coordinate space. 



3.2 Chebyshev basis representation 

The long term dynamics are thus contained in the operator (J — B)^ 1 and standard properties of Perron- 
Frobenius operators ensure that the sum over n in equation (Q converges for non-vanishing dissipation. 
In order to evaluate (/ — B)~ l it is convenient to express the operator B in a suitable set of basis 
functions defined on the boundary. In Ref. [7] a Fourier basis has been applied, which is a natural choice 
of a complete basis for problems with periodic boundary conditions. However, a number of difficulties 
arise with this choice such as slow convergence of quadrature rules for the associated integrals and the 
treatment of corners on the boundary. 

Here we employ a Chebyshev polynomial basis representation with Gauss- Chebyshev quadrature, in 
which case the integration is optimal for polynomial-type integrands. Problems due to singular behavior 
at corners are avoided due to integrating over phase space, rather than over pairs of boundary coordinates. 
Gauss- Chebyshev quadrature incorporates the orthogonality weight functions for the Chebyshev basis 
automatically in the quadrature rule. In the case d = 2, we have s G [0, L) and Chebyshev basis may be 
expressed in the form 

Tn(X S ) = X [^-T ni ( 2 4 - lW 2 feV (10) 



\p\L ™\L J n2 \\p\ 

with n = (711,712) non-negative integers and T ni the Chebyshev polynomial of order n\. The Chebyshev 
basis approximation B mn of B may be written: 

r W m (X s )f m (X s )B(X s , Y s - u)f n (Y s )dY s dX s , , 

W m {r{Y s ))f m {r{Y s ))w{Y s )f n {Y s )dY\ 



where <9P = [0, L) x (— \p\, \p\) is the phase space on the boundary at fixed "energy" H{X) = \p\ 2 and 

W (X s ) = ^ mi ^ m2 1 1 (12) 

* 2 a/i - wm - 1) 2 v / i r w' 

is the weight function for the inner product in which the Chebyshev basis is orthonormal. Here 70 = 1/2 
and 7 ni = 1 for n\ — 1, 2, .... It is convenient for the Gauss-Chebyshev quadrature if the argument in the 
weight function W m is the same as the integration variable and so a change of variables Y s = ^ W (X S ) 
with ip w = (0 W ) _1 is carried out to give 

w m {x s )f m {x s ) w {r{x s ))f n {r{x s ))\j{x s )\dx s . (13) 
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Here the Jacobian term is 



|J(X*)|=|^(X*)| 



d£_ 

ds 

ds 



ds' 

dps 

dp's 



(14) 



which is equal to one for Hamiltonian flows and takes account of changes in the wavenumber between 
subsystems. In this representation the integration is with respect to position and momentum at the end 
of the ray being considered, and we are mapping back to the initial point using ip u . 



3.3 Subsystems 

Recall the splitting into subsystems i = 1, .., Nq introduced earlier. The dynamics in each subsystem 
are considered separately so that both variability in the wave velocity q and non-convex domains may 
be handled simply. Coupling between sub-elements can then be treated as losses in one subsystem and 
source terms in another. Typical subsystem interfaces are surfaces of reflection/ transmission due to 
sudden changes in material parameters or local boundary conditions. We describe the full dynamics in 
terms of subsystem boundary operators By] flow between Qj and fli is only possible if fij fl flj ^ and 
one obtains 

%(x*,x;) = Wij (xs)6(x* - 4%{xf)), (is) 

where <f>fj is the boundary map in subsystem j mapped onto the boundary of the adjacent subsystem 
i and Xf are the boundary coordinates of The weight contains, among other factors, reflection 
and transmission coefficients characterizing the coupling at the interface between Qj and 

A basis function representation BJp of the full operator B as suggested in equation (TIT]) is now 
written in terms of subsystem boundary basis functions T % n to give 

B$ n =[ [ W m (Xt)fl{Xt)B{Xt,X^u)f^)dX^dXt. (16) 

Here <9Pj is simply the boundary coordinate phase space for The equilibrium distribution on the 
interfaces of the subsystems is then obtained by solving the system of equations (JS} 

(I-B)p r = pP. (17) 

Here B is the full operator including all subsystems and the equation is solved for the unknown energy 
densities pr = (prji=i,..,7v n , where pr, denotes the (Chebyshev) coefficients of the density on Tj, the 
boundary of fV 



4 NUMERICAL IMPLEMENTATION 

4.1 From SEA to DEA 

Up to now, the various representations given are all equivalent and correspond to a description of the wave 
dynamics in terms of the ray tracing ansatz (JSJ). Traditional ray tracing based on sampling ray solutions 
over the available phase space is rather inefficient. Convergence tends to be fairly slow, especially if the 
absorption is low and an exponentially increasing number of long paths including multiple reflections 
need to be taken into account. 

An SEA treatment emerges when approximating the individual operators By in terms of constant 
functions only [7]; using, for example, a Fourier basis this would correspond to an approximation in 
terms of the lowest order basis functions only. In the case of a Chebyshev basis, the integrand in (flB]) 
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Figure 1: (Color online) A polygonal configuration with Tj = dQj for j = 1,2 showing further subdivision 
Tj. of the boundary. The interface is formed by the boundary sections = T 2l . A separate spatial 
basis function approximation is applied for each subdivision. 

is not constant due to the presence of the weight function W m and an SEA treatment is obtained only 
after restricting the basis to Tq and omitting the weights Wq. Note that when using Gauss- Chebyshev 
quadrature an explicit division by Wo is required since the weight function is automatically included in 
the quadrature rule. In the SEA case the matrix By is one-dimensional and gives the mean transmission 
rate from subsystem j to subsystem i. It is thus equivalent to the coupling loss factor used in standard 
SEA equations [2]. The resulting full TV^-dimensional B matrix yields a set of SEA equations using the 
relation ([17]) after mapping the boundary densities back into the interior. 

An SEA approximation is justified if the ray dynamics within each subsystem are sufficiently chaotic 
such that a trajectory entering subsystem j forgets everything about its past before exiting Qj] SEA can 
thus be viewed as a Markov approximation of the deterministic dynamics. Thus correlations within the 
dynamics must decay rapidly on the time scale it takes for a typical ray to leave Qj. This condition will 
often be fulfilled if the subsystem boundaries are sufficiently irregular, the subsystems are dynamically 
well separated and absorption and dissipation is small, conditions typically cited in an SEA context. In 
this case SEA is an extremely efficient method compared to standard ray tracing. However for subsystems 
with regular features, such as rectangular cavities or corridor-like elements, incoming rays are directly 
channeled into outgoing rays, thus rendering the Markov approximation invalid and introducing memory 
effects. Likewise, strong damping may lead to a significant decay of the signal before reaching the exit 
channel introducing geometric (system dependent) effects. 

The features that cause an SEA approximation to fail as described above are incorporated into 
the model here by including higher order basis functions and weight functions W m for each subsystem 
boundary operator By. It therefore becomes possible to resolve the fine structure of the dynamics and 
their correlation along with effects due to non-uniform damping over typical scales of the subsystem. 
The maximal number of basis functions needed to reach convergence is expected to be relatively small 
thus making the new method more efficient than a full ray tracing treatment, particularly when damping 
is low. 

Representing the ray dynamics in terms of finite dimensional transition matrices corresponds to a 
refinement of an SEA technique, taking advantage of the efficiency of SEA but including information 
about the ray dynamics when necessary. It overcomes some of the limitations of SEA and puts the 
underlying SEA assumptions on sound foundations. 

4.2 Chebyshev DEA 

A favourable property of the Chebyshev basis approximation compared with a Fourier basis is the flexi- 
bility it allows in the choice of approximation regions due to not requiring periodic boundary conditions 
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for convergence. In many cases it will be advantageous to subdivide the subsystem boundaries and apply 
a Chebyshev basis for each subdivision. An example where this subdivision takes place at the vertices 
of each subsystem is shown in Fig. [TJ A separate spatial basis expansion is then applied in a number 
of sections Fj v i — 1, ...,Nj of the boundary Tj of a particular subsystem Qj, j = 1, ..,Nq. This leads 
to separate entries in the matrix representation B for each boundary section Tj. and hence a larger but 
sparser matrix. 

In order to understand why this geometric subdivision may be beneficial, we first need to consider 
some properties of the basis function expansions. The basis expansion of the boundary operator B 
results in density functions and pr t in equation fl9]) being expanded in terms of the chosen basis 
functions. The convergence rate of the basis expansion therefore depends on the properties of these 
density functions and in particular their smoothness [15] . The basis function approximations are carried 
out with respect to both position s and momentum along the boundary p s at the endpoint of the ray, 
recall equation (fT5|) . A discontinuity in the normal to Tj for some i — 1, ...,Nq will, in general, result 
in a non-differentiable initial density. It is therefore recommended to employ separate Chebyshev basis 
approximations with respect to s for regions of T split by an edge or vertex, see Fig. [TJ It is also 
recommended to employ a separate approximation with respect to s where there is a sudden change in 
boundary conditions resulting in a non-smooth initial density, for example, along an interface between 
two subsystems. The geometric subdivision is therefore to enable smoothness of the density functions, 
Py^ and pr\ , over the regions of approximation and therefore (geometric) exponential order convergence 
in the Chebyshev basis expansions. 

The method described above has been implemented numerically for a variety of problems with Q C 
R 2 and homogeneous Dirichlet boundary conditions for the total wave, (that is, including the initial 

— 1/2 

contribution from the source) on the outer boundary T. For simplicity we consider = g { for all 
i = l,..,Nci, that is, we set the adiabatic compressibility to unity. Results are compared with the 
numerically exact solutions obtained from state-of-the-art adaptive FEM software. A brief account of 
the FE method is given in the next section. 

Let s G [0, LjJ and p s G (—ki, parameterize the associated phase space where Lj t is the length of 
Tj i . In these local sub-system boundary coordinates the Chebyshev basis is given by 



<«-/U J^f-T ni (s)T n2 (p s ), (18) 

2 Ji 

where s = (2s /Lj.) — 1 and p s = p s /k{. The weight function W n is given by 

W n (s, Ps ) = 4 ^ * -=L=. (19) 
vr 2 VI - s 2 a/1 - p 2 s 

The transmission probability at the intersection of two subsystems is given by [20] 

/ t i n \ ^{kout/kiYi ) cos(6» in ) cos(6» out ) 

( (kout/ km ) cos(# in ) + cos(6Ui)) 2 ' 

with 6 in {6 ou t) the angle between the incoming (outgoing) ray and the inward normal to the boundary 
and k in (k out ) the wavenumber in the subsystem through which the incoming (outgoing) ray is traveling. 
Incoming and outgoing rays are related through Snell's law {k in sm(9 in ) = k out sm(9 out )), and hence 9 0Ut 
may be calculated from the other three quantities. 

Given the end point s G Tj t and the incoming momentum p s , one can obtain the corresponding ray. 
With sufficient geometric knowledge of Tj and the boundaries of its adjoining subsystems, the start point 
of the ray s' can be determined as its intersection with one of these boundaries, say in the subset Tp a 
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for some j3 a = 1, N a and a = 1, Nq. Once this is known it is straightforward to obtain the initial 
momentum p' s . Writing out the Jacobian ( 1T4|) one obtains 

B nl = T I I W m (s,p s )fi(s,p s )f^(s',p' s )x 

WjiP a ( s ',P' s )dsdp s , 

where Wj^ a (s' ,p' s ) = exp(—fiiL)wJ ij3ct (s',p' s ). Here /ij is the damping coefficient in Qi as before, L is the 
length of the trajectory from s' to s and the reflection/ transmission coefficients are 

r / / / \ J 3ia if s ^ Ti a 

WjA{S > p ' } -\ S ia + (-lY^wtih, k a , 6$ if s e r ia . 

Here, T ia denotes any subset of T a forming the intersection of two subsystems (that is, for i ^ a, 
T ia = Ti n r a ) and S ia is the Kronecker delta. Also k™ in = ^sin(^™ n ) and k™ ax = ^ sin(#™ a:c ) are 
the minimum and maximum values of p s , respectively, where 6™ m , 6™ ax G (— 7r/2,7r/2) are the angles 
between the inward normal to Tj t and the rays from s' to each of the ends of Tj i . 

All DEA/ SEA computations are performed on a desktop PC with a dual core 2.83 GHz processor 
using C++ with Diffpack (www.diffpack.com). 

4.3 /ip-adaptive discontinuous Galerkin finite element method 

In the Sec. refseexomputations, we will compare DEA results twith numerically exact solutions of the 
wave equation ([T]). Discontinuous Galerkin (DG) methods have become increasingly popular for elliptic 
problems in recent years [21] and is also our method of choice for calculating the Green function. The 
main reason for this interest in DG methods is that allowing for discontinuities across elements gives 
extraordinary flexibility in terms of mesh design and choice of shape functions. Additionally, /ip-adaptive 
DG methods, which are based on locally refined meshes and variable approximation orders, achieve 
tremendous gains in computational efficiency for challenging problems [221 ESI [2H 125] ■ 
In the following, we consider the wave equation 

(c^A + c^)G(r, r ; u>) = -8(r - r ), i = 1, N n , (22) 

with uji = w+z/ijCj/2 and refljC M 2 . The Green function G is related to G in equation (JTJ by G = c^G, 
where Co is the wave velocity in the subsystem Q$ containing the source point r$. 

To drive the hp- adapt ivity we use an explicit energy norm a posteriori error estimator inspired by 
Refs. [26| [27] . We apply a simple fixed- fraction strategy on the error estimator to mark the elements 
to adapt. For each marked element, the choice of whether to locally refine it or vary its approximation 
order is made by estimating the decay of the coefficients in an L 2 -orthogonal polynomial expansion in 
order to test the local analyticity of the solution in the interior of the element |28j. The details of the 
implementation of the method can be found in Appendix |A] 

All finite element simulations have been carried out using AptoFEM (www . aptof em . com) on a parallel 
machine. 

5 COMPUTATIONAL RESULTS 

5.1 Coupled two-cavity systems 

A variety of two-cavity systems are considered as in Ref. [7] and are shown in Fig. [2J The coordinates of 
the vertices and source points of these systems are given in Appendix [Bl We compare our DEA results 
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Figure 2: Coupled two-domain systems: configurations A, B and C, respectively. 
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Figure 3: (Color online) Ratio of total energies R = ||Gi|| 2 /||G2|| 2 m configuration A with c x = c 2 = 1 
(left) and c\ = 0.5, c 2 = 1 (right). Dotted lines show FEM results computed within ±5Hz of the DEA 
results. 

with standard SEA calculations; relevant parameters such as the area of each subsystem have been listed 
in the Appendix [Bl The SEA coupling loss factors are obtained from the transmission coefficients (1201) 
integrated over the wavenumber ki n of the incoming rays. The attenuation is given by the damping factor 
//; in the following we will assume hysteretic damping with /ij = ur// (2q) for i = 1,2. 

Configuration A features irregular shaped well separated pentagonal subsystems and thus SEA is 
expected to work well. In configuration B the size of the interface between the subsystems is increased 
reducing their dynamical separation and therefore the applicability of SEA. Configuration C includes a 
rectangular left-hand subsystem channeling rays out of the subsystem and introducing long-range corre- 
lations in the dynamics. In addition, the source is further from the intersection of the two subsystems. 
SEA is thus not expected to work well for this configuration. Note that SEA results are in general 
insensitive to the position of the source, whereas actual trajectory calculations may well depend on the 
exact position. 

Finite basis sets have been employed with ni,n,2 = 0,..,N, which gives rise to matrices of size 
dimi? = (N + l) 2 (iVi + AT 2 ), with basis functions of the same order for position and momentum in both 
subsystems. Note that Ni,i = 1,2 denotes the number of subdivisions of the boundary of subsystem i 
as defined in Section 14.21 Energy distributions have been studied as a function of the frequency with a 
hysteretic damping factor rj = 0.01. Here and in the remainder of this work the subsystems are numbered 
1, Nq from left to right. 

Fig. [3] shows approximations of the ratios of the total energy in each subsystem ||G r i|| 2 /||G , 2|| 2 in 
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Figure 4: (Color online) Ratio of total energies R = ||Gi|| 2 /||G2|| 2 in configuration B (left) and configu- 
ration C (right). Dotted lines show FEM results computed within ±5Hz of the DEA results. 



configuration A where 

||G 2 ,|| 2 := / \G{r,r Q] u)\ 2 dr, i = l,...,N Q . (23) 

The left subplot shows the results with c\ = = lms -1 and the right with c\ = 0.5ms -1 and c<i = 1ms -1 . 
The dotted lines represent solutions computed using the FEM (described in section |4"73"|) at an equi-spaced 
range of frequencies within ±5Hz of the frequencies used for the SEA and DEA computations. Note that 
the damping is fixed to the value employed for the central (SEA/DEA) frequency. In the right-hand 
plot the results only go up to 50Hz due to the high computational cost of the FEM code for the large 
wavenumbers in subsystem 1. Fast convergence with increasing basis size is evident in each case. It is 
clear that SEA works reasonably well for configuration A since the SEA prediction is close to that from 
DEA and within the range of FEM solution values. In particular the SEA prediction is good for low 
damping values (that is, low frequencies). At closer inspection one notices also that the divergence from 
the DEA result is smaller in the case when c\ ^ C2, thus demonstrating an increase in the dynamical 
separation between subsystems in this case. Comparing SEA with the N = Q case between 10Hz and 
50Hz, the SEA results differ from DEA by between 8% and 15% when c\ = c<i = lms -1 , but only by 
between 5% and 7% when c\ = 0.5ms -1 and C2 = lms -1 . 

Figure H] shows approximations of the ratios of the total energy in each subsystem in configurations 
B and C with c\ = C2 = lms -1 . The dotted lines are as before and fast convergence with increasing 
basis size is again evident, although a slightly higher order was required for configuration C to capture 
the exponential decay due to dissipation along the rectangular cavity. As expected, the SEA prediction 
diverges from both the FEM and DEA predictions in configurations B and C. Comparing SEA with the 
N = 6 case for configuration B between 10Hz and 50Hz (for consistency with the data for configuration 
A), the SEA results differ from DEA by between 18% and 29%. For configuration C the deviation is 
about 30% in the range 10Hz to 50 Hz, although this grows considerably larger if one includes the data 
for 60 and 70Hz. 
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Figure 5: (Color online) Distribution of log 10 (|G| 2 ) in a five cavity system. 



5.2 Complex built-up systems 

In this section the versatility and efficiency of the Chebyshev approximation with Gauss-Chebyshev 
quadrature are demonstrated by considering a complex built-up system. A configuration of five coupled 
acoustic cavities is considered as shown in Figure El The coordinates of the vertices and source point 
are given in Appendix [Bj The solution in the interior of each cavity is plotted and the source point 
in the central cavity is clearly evident. The subsystems are each convex polygonal regions as before, 
the jagged appearance of the boundary is a result of the plotted region being formed of the largest 
regular grid fitting strictly inside the boundary of each subsystem. The wave velocities are taken to 
be c\ = C2 = C4 = C5 = lms -1 and C3 = 0.5ms _1 . Figure |5] shows the DEA approximation of the 
distribution of log 10 (|G| 2 ) throughout the system with N = 8. The logarithm of the solution is taken 
since a large range of values are present between the peak of the source and the extremal subsystems. 
The plot is for the 10Hz case with the same frequency and damping correspondence as in the earlier 
two-cavity configurations. DEA clearly gives much more detailed spatial information about the wave 
energy distributions than SEA, which assumes a constant density in each subsystem. In particular here 
one can see local variations close to subsystem interfaces and a drop in the intensity as one moves away 
from the source. Note that more energy flows into the far right subsystem as compared to the far left 
subsystem due to there being a direct channel for the energy to travel along to the right of the source, 
but not to the left. 

Figure E] shows approximations of ||Gj|| 2 for i = 1,..,5 computed using both SEA and DEA up to 
N = 8. The three figures represent the total energy (actually Cq = 1/16 multiplied by this quantity 
since G rather than G has been computed) obtained for each subregion at three different frequencies 
and thus damping levels; using the same parameters as before the 10, 20 and 30 Hz cases are shown 
from left to right. Note that results are shown on a logarithmic scale and that the overall amplitude 
decreases with increasing frequency due to increased damping. Here a comparison with FEM simulations 
is not considered due to the high computational cost associated with computing high and multi-frequency 
solutions over large domains. We can deduce that SEA is working very well in subsystems 1 to 3 due 
to the level of agreement with the DEA calculations. Very few trajectories can move from the source 
directly into subsystem 1 and multiple scattering events lead to a local equidistribution, that is, incoming 
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Figure 6: (Color online) A plot of ||Gj|| 2 for subsystems % = 1, ..,5 at three different frequencies. 



and outgoing rays in subsystem 2 are uncorrelated. The situation is different for subsystems 4 and 5 
where the influence of the direct channel from the source to subsystem 5 becomes important. Compared 
to the SEA result, the DEA calculations show a noticeable increase in the values of 1 1 1 1 2 and a slight 
decrease in the values of HG4II 2 , that is, more energy reaches subsystem 5 than predicted by a Markov 
approximation of the dynamics. The DEA calculations again appear to converge reasonably quickly 
between N = 6 and iV = 8. It is also evident that SEA works best for the lower damping values here, 
which is consistent with our observations for the two-cavity configurations. 



6 CONCLUSIONS 

Dynamical energy analysis has been reformulated in terms of a Chebyshev basis expansion allowing a 
greater degree of flexibility in the design of the method and more efficient coding. A comparison of the 
numerical results with finite element method computations showed that dynamical energy analysis is 
robust and retains accuracy in cases when SEA fails. Examples where the wave velocity varies between 
subsystems were also considered for the first time using DEA. An extension to built-up systems, facilitated 
by the efficient Chebyshev basis reformulation of DEA, showed that DEA is a flexible, robust and efficient 
method for estimating energy density distributions at high frequencies. 



Acknowledgments 

Support from the EPSRC (grant EP/F069391/1) and the EU (FP7IAPP grant MIDEA) is gratefully 
acknowledged. 

The authors also wish to thank Inutech Gmbh, Niirnberg for Diffpack guidance and licences, and 
Dmitrii Maksimov for carefully reading the manuscript. 



A Implementation of the DG method 

In order to compute an approximation of G in (1221) . we partition the domain Q with shape- regular 
meshes Th formed by open triangles {K}k£T h - We assume that in the interior of each element K e Th, 
the wave velocity and the damping are constant, such values are denoted in the interior of K by ck 
and hk- Further, each element K can then be affinely mapped onto the generic reference element K. 
The diameter of an element K e Th is denoted by hx- Due to our assumptions on the meshes, these 
diameters are of bounded variation, that is, there is a constant b\ > 1 such that 

< h K /h K , < 61, (24) 
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whenever K and K' share a common edge. We store the elemental diameters in the mesh size vector 
h = {hx '■ K G %,}■ Similarly, we associate with each element K G Th a polynomial degree px > 1 and 
define the degree vector p = {px '■ K <E %,}■ We assume that p is of bounded variation as well, that is, 
there is a constant b 2 > 1 such that 

62 1 < VkIvk< < 62, (25) 

whenever X and iT' share a common edge. 

For a partition 7/i of f2 and a degree vector p, we define the hp-version discontinuous Galerkin finite 
element space Vh of complex valued functions by 

V h = {veL 2 (Q) : v \ K eV PK (K), KeT h }. (26) 

Here, V PK (K) is the space of polynomials on K of total degree less than or equal to pk- 

Next, we define some trace operators that are required for the DG methods. To this end, we denote 
by £x(7h) the set of all interior edges and by £r{Th) the set of all boundary edges of the partition Th- 
Furthermore, we define £(%,) = £i(lh) U£r(7/i). The boundary OK of an element K and the sets dK\T 
and dK D T will be identified in a natural way with the corresponding subsets of £(%,)■ 

Let K + and K~ be two adjacent elements of 7^, and k G £x(Th) given by k = dK + fl dK~ . 
Furthermore, let v be a complex scalar- valued function, that is smooth inside each element K ± . By v ± , 
we denote the traces of v on k taken from within the interior of K , respectively. Then, the weighted 
average of the diffusive flux c 2 VhV along k G £x(7h) is given by 

ic 2 V h v}} = c K-^f + cl + c 2 K -V h v^ 

C K+ + C K~ 

Similarly, the jump of v across k G £xiTh) is given by 

\v\ = v + n K + + v~ n K - , 

where we denote by n K ± the unit outward normal vector of dK^, respectively. 

On a boundary edge k G £r{Th), we set {[c^V/iv]} = c 2 VhV and [v] = vn , with n denoting the unit 
outward normal vector on the boundary Y. 

For a mesh Th on fl and a polynomial degree vector p, let Vh be the hp-version finite element space 
defined in (|26l) . We consider the (symmetric) weighted interior penalty discretization [29] of ( 122]) : find 
(5^ G V/t such that 

= F fc («) , for all v G V h , (27) 



where 



/2) wdr 

- E / (I^W • [ u ] + {c 2 V h u} ■ Jv\) ds 

-Ph(f) := / #(r - r )vdr , 

and Vfe denotes the element wise gradient operator. Since DG methods allow for discontinuities in the 
finite element approximation, the supports of the shape functions never extend to more than one element. 
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A consequence is that the stencil is minimal in the sense that each element communicates only with its 
direct neighbors and the communication happens across the edges of the mesh where the functions in 
the finite element space are not continuous. This makes it natural to include the last two terms in 
the definition of Ah(-,-), which control the average and the flux of discontinuous functions across the 
edges. In particular the stability of the method is guaranteed by the third term, which penalizes the 
discontinuities across the edges. So the second and the third terms are peculiar to the DG method of 
choice, in contrast to the first term which depends on the PDE under consideration. Furthermore, the 
function c G L°°(£(Th)) is the discontinuity stabilization function that is chosen as follows: we define 
the functions h G L°° (£(%)) and p G L°° (£(%)) by 



h(r) : = 
p(r) := 



min(h K , h K >), r G k G £x{Th), k = dK D dK' , 

h K , r G k G £r(7/i), k> G dK D T, 

m&x(p K ,p K >), r G k G £x(Th), k = dK n dK', 

[p K , r G k G £ r (T h ), k G dK n T, 



and set 



r 2 r 2 d 2 

„, C K C K' P / OQ \ 

7 c 2 +c 2 "h' (28) 

with a parameter 7 > that is independent of h, p, ck and ck< ■ The definition of c is equivalent 
to an hp- version of the weighted penalty parameter in Ref. [29]. The definitions of h(r), p(r) and c 
are designed to work with /ip-adaptivity where adjacent elements of very different sizes and orders are 
possible. In order to keep the method stable under the adaptation process it is necessary to penalize 
more discontinuities across the smaller edges and also across higher order elements. This can be seen 
straightaway from 



B Geometric and SEA data 

For completeness the coordinates of the vertices for the systems treated numerically in the paper are 
detailed in Table [H For consistency with the units quoted in the paper the distances should be considered 
in metres. In each configuration the vertices are ordered anticlockwise starting from the vertex with the 
maximum x-coordinate. For configurations A and B the source point is (—0.4,0.5), for configuration 
C it is (—1.4,0.4993) and for the five plate configuration the source is located at (—0.5,0.4993). In 
configurations A, B and C the interface between the subsystems is located on the intersection with the 
line x = 0. For the five cavity configuration the interfaces are on the intersections with the lines x = —2, 
x = —1, x = and x = 1. Table Ogives additional geometric parameters useful for SEA verification. 
Note that in a DEA framework the SEA coupling loss factors are computed using the transmission law 
(120]) and details of the damping, frequency and wave velocity are provided in Section HI 
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